Long-range-corrected hybrids including RPA correlation 
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We recently demonstrated a connection between the random phase approximation (RPA) and 
coupled cluster theory [J. Chem. Phys. 129, 231101 (2008)]. Based on this result, we here propose 
and test a simple scheme for introducing long-range RPA correlation into density functional theory. 
Our method provides good thermochemical results and models van der Waals interactions accurately. 



We recently demonstrated that the ground state corre- 
lation energy associated with the random phase approx- 
imation (RPA)L2.2 i M i gi l £ i s connected with an approxi- 
mate ring coupled cluster doubles (rCCD) approach. The 
RPA excitation problem requires the solution of 
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where the matrices A, B, X, and Y are of dimension 
ov x ov, with o and v being the number of occupied and 
unoccupied spin-orbitals, respectively. The plasmonic 
formula^ for the RPA ground state correlation energy 
is 
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RPA 
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-Tr(w- A). 
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As shown in our previous paper ^ Eq. [T]is equivalent to 



B + AT + TA + TBT = 0, 



(3) 



an approximate coupled cluster (CC) doubles equation 
with excitation amplitudes T = YX" 1 . Further, the 
RPA correlation energy can be evaluated from the CC- 
like expression 
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-Tr(B T) = -Tr(w - A 
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The excitation amplitudes T imply the existence of 
an underlying wave function in RPA. This lets us fol- 
low Savin, Stoll, and coworker s 9 ' 10 and introduce range 
separation between density functional theory (DFT) 
and wave function theory, with RPA for the long- 
range correlation. Range separation is a powerful 
tcchniqu o 11 ' 12 ' 13 ' 14 ' 15 ' 16 ' 17 ' 18 ' 19 that can improve upon 
both standard wave function methods and semilocal 
exchange-correlation (xc) functionals. Range separated 
DFT partitions the electron-electron interaction opera- 
tor into short (SR) and long (LR) ranges 



fc(wri 2 ) erf(wri2) 
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and typically (but not always, see Refs. Il3lfl9h treats 
the SR (LR) component with semilocal (wave function) 
approximations. 



Toulouse et aZ.— recently proposed a range-separated 
treatment of RPA via the adiabatic-connection 
fluctuation-dissipation theorem,— combining a 
short-range semilocal xc functional with long-range 
full RPA. Consistent with other work ; 17 ' 18 they find 
that their long-range correlation energy has a relatively 
weak basis set dependence. They also remove some 
artifacts of full-range RPA, including a "bump" in the 
symmetry-restricted singlet Bc2 dissociation curve. 

Theory. We propose a simpler long-range RPA based 
on the connection to coupled cluster theory discussed 
above, which avoids a costly adiabatic connection inte- 
gral. We evaluate long-range RPA as a one-shot correc- 
tion to a self-consistent generalized Kohn-Sham^ (GKS) 
calculation combining long-range exact (Hartree-Fock- 
type, HF) exchange and short-range local spin density 
(LSDA) exchange-correlation. Our xc energy is 
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We evaluate E^ R ~ RPA with Eq. [30 and build A and B 
from the long-range two-electron integrals and the GKS 
spin-orbitals and orbital energies. The coefficient crpa 
is discussed below. 

In this work, we focus on what we will refer to as direct 
RPA. In the (real) canonical spin-orbital basis we use 
throughout this paper, the direct RPA matrices are 

A ia ,jb = {e a - ei)6ij5 ab + (ib\ aj) , (7a) 
Bia,jb = ab) . (7b) 

Here e is a generalized Kohn-Sham orbital energy In- 
dices i and j indicate occupied spin-orbitals, a and b in- 
dicate virtual spin-orbitals, and ( ij \ ab) is a two-electron 
integral in Dirac's notation. For real orbitals, note that 
( ib | aj) = (ij\ ab). What we refer to as full RPA uses 
antisymmetrized two-electron integrals in Eq. [7J 

For practical calculations, direct RPA has the great 
advantage that the correlation energy is guaranteed to 
be real if the orbitals obey the aufbau principle. This 
is not true of full RPA, where instabilities in the refer- 
ence determinant^ can produce a complex correlation 
energy^ In cases such as Be2 where the reference has 
triplet instabilities^ one may have to limit full RPA to 
singlet excitations, while with direct RPA no such re- 
striction is needed. 

An additional advantage of direct RPA is that it re- 
duces the dimension of the problem compared to full RPA 
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(ring CCD with antisymmctrizcd two-electron integrals). 
Consider the block of B corresponding to spatial orbitals 
(pi, ipA; <fj, and ifB and spin ordering aa, (3(3, a/3, and 
(3a. The B matrix for full RPA and closed shells becomes 



N atom total energies 
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with J = (IJ\ AB) and K = {IJ\ BA) . For open shells 
(unrestricted), the J and K, entries in different blocks will 
generally differ. Direct RPA (K, = 0) zeros the "spin-flip" 
block of B containing only —K, reducing the dimension 
by a factor of 2 compared to full RPA. For both closed 
and open shell systems, we can diagonalizc the remaining 
upper block of B with the unitary transformation 
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For closed shells, Q is a spin-adaptation operator. The 
resulting 2x2 diagonal matrix for direct RPA has only 
one non-zero eigenvalue, corresponding to singlet excita- 
tions. Since only the singlet block of B is non-zero, we 
only need the singlet part of T to evaluate the correlation 
energy of Eq. [4] Additionally, in blocks where B = the 
ring CCD equation (Eq. [3]) becomes 



AT + TA = 0. 



(10) 



Because the direct RPA A is positive definite,— the only 
solution is T = 0. Thus the spin flip components of T 
vanish, and in the closed shell case all triplet components 
of T vanish. Triplets do not contribute to direct RPA. 

A third advantage of direct RPA is that the correla- 
tion energy expression is well defined, unlike in full RPA. 1 
While in both direct and full RPA, wc have Tr(B T) = 
Tr(u; — A),— only in direct RPA is the prcfactor in the 
correlation energy unambiguously 1/2. For full RPA, the 
plasmonic formula suggests the prefactor should be 1/2, 
while the connection to CC theory suggests that the pref- 
actor should be 1/4. Jcziorski and coworkers have pre- 
sented another alternative full RPA energy expression 
that uses the plasmonic prefactor of 1/2, but subtracts 
the second-order MP2 correlation energy^ 

The elimination of exchange integrals in direct RPA 
may appear artificial from a wave function perspective, 
and in fact the wave function underlying direct RPA can 
violate the Pauli exclusion principle. However, direct 
RPA can be rigorously derived by applying the adia- 
batic connection fluctuation-dissipation theorem to the 
Kohn-Sham noninteracting reference system4£ On bal- 
ance, long-range direct RPA from the GKS reference ap- 
pears to be a practical ansatz for long-range correlation. 

Computational details. We have implemented our ex- 
pressions into the development version of the GAUSSIAN 
suite of programs^ Matrices A and B are evaluated in 
the full spin-orbital basis set, without symmetry adap- 
tation. We use the range-separated LSDA correlation 
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FIG. 1: Basis set dependence of LC-cjLDA+dRPA. (Top) To- 
tal energy of N atom. (Bottom) Dissociation energy of N2 at 
experimental bond length 1.098 Angstrom. Results are calcu- 
lated with aug-cc-pV^Z basis sets, and normalized to the aug- 
cc-pV5Z result. LC-tjLDA+dRPA calculations use chpa — 1 
and either uj = 0.5 or lu = 1.2 Bohr -1 . 



functional of Paziani and coworkers.— In what follows, 
"LC-cjLDA" denotes GKS calculations using short-range 
LSDA xc and long-range HF exchange. (Note that in 
our previous papers, this acronym denoted SR LSDA 
exchange, LR HF exchange, and full range LSDA cor- 
relation.) Adding long-range direct RPA correlation 
as described above results in "LC-wLDA+dRPA" . The 
"dRPA" acronym by itself denotes conventional, full 
range (i.e., not range separated) HF exchange and di- 
rect RPA correlation; this "dRPA" energy is evaluated 
from self-consistent Kohn-Sham orbitals and orbital ener- 
gies calculated with the Perdew-Burkc-Ernzerhof (PBE) 
generalized gradient (GGA) xc functional.— Open-shell 
systems are treated spin unrestricted. Eq. [3] is solved it- 
eratively using DBS 29 ' 30 for coupled clusters. Correlated 
calculations use frozen core electrons. 

Like other workers ^^sa^A^x^vq we select 
the range separation parameter 10 empirically. Our stan- 
dard LC-wLDA+dRPA calculations set crpa = 1 in Eq. 
[6] We also explore treating crpa as an empirical param- 
eter. 
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TABLE I: Mean absolute errors (kcal/mol) in AE6 at- 
omization energies, G2/97 heats of formation, BH6 and 
HTBH38/04 hydrogen-transfer reaction barrier heights, and 
NHTBH38/04 non-hydrogen-transfer barrier heights. 6- 
311+G(2d,2p) basis set. lo in Bohr -1 . 
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FIG. 2: Counterpoise-corrected dissociation curves of van der 
Waals dimers. Aug-cc-pV5Z basis, other details in the text. 

Numerical results. Figure [1] shows the basis set de- 
pendence of LC-wLDA+dRPA for the total energy of N 
atom (top) and the dissociation energy of N2 (bottom) . 
As in previous work ) 17 i 20 i 32 the long-range correlation has 
a much weaker basis set dependence than full-range cor- 
relation. Calculations at the relatively large lo = 1.2 
Bohr -1 , which incorporate a larger fraction of direct RPA 
correlation, converge more slowly with basis set size. 

Figure [2] shows counterpoise-corrected dissociation 
curves of van der Waals dimers He2, Ne2, and Ar2, 
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FIG. 3: Dissociation of symmetry-restricted singlet H2. Aug- 
cc-pVTZ basis set, lo = 1.2 Bohr -1 . "Accurate" results are 
full CI/aug-cc-pV5Z. 



evaluated in the large aug-cc-pV5Z basis setM- Accu- 
rate curves arc from Ref. [3J. LC-wLDA and LC- 
wLDA+dRPA use lo = 1.2 Bohr -1 . Rescaling the long- 
range dRPA correlation with crpa = 1.5 significantly im- 
proves the results, suggesting that the rescaling primarily 
corrects for beyond-dRPA correlation effects rather than 
basis set incompleteness. The weak He2 binding of stan- 
dard (full-range) dRPA is increased by using BP86 rather 
than PBE orbitals and orbital energies, consistent with 
Ref. (not shown). 

Table U shows mean absolute errors in the small AE6 
and BH6 sets of 6 atomization energies and 6 reac- 
tion barrier heights, the G2/97 set of 148 heats of 
formation 3 ^ and the HTBH38/04 and NHTBH38/04 sets 
of 38 hydrogen-transfer and 38 non-hydrogen-transfcr 
barrier heights . 37 ' 38 Calculations use the 6-311+G(2d,2p) 
basis set. G2/97 calculations use B3LYP/6-31G(2df,p) 
geometries and vibrational frequencies^ other geome- 
tries and reference values are taken from Refs. I35ll36l . 
I37ll38l . Results are presented for the thermochemically 
optimal lo. LC-wLDA+dRPA is quite accurate for ther- 
mochemistry and kinetics, particularly with the cmpir- 
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ical rescaling c^pa = 1.5. It improves upon the un- 
derlying LC-wLDA, giving results comparable to the ac- 
curate LC-ojPBE combination of range-separated GGA 
exchange and full-range GGA correlation^ Adding long- 
range dRPA to LC-wPBE appears to give significant 
double-counting of correlation (not shown). However, 
combining long-range dRPA with a short-range GGA (as 
in Ref. 1201 ) may provide further improvements. 

One of our goals is to apply range-separated dRPA 
to metallic systems. This is possible because dRPA is 
robust as the band gap closes. Figure [3] illustrates LC- 
wLDA+dRPA in one such prototypical system. The fig- 
ure plots the energy of spin- and symmetry-restricted 
aug-cc-pVTZ H2 as a function of H-H bond length. The 
MP2 correlation energy diverges as the bond length in- 
creases and the HOMO-LUMO gap approaches zero. Full 
RPA (not shown) yields non-real correlation energies for 
bonds stretched beyond the Coulson-Fischer point. In 
contrast, the LC-cjLDA+dRPA energy is real and finite 
at all bond lengths. Long-range dRPA does not capture 



all of the nondynamical correlation present in stretched 
symmetric singlet H 2 , and the dissociation limit is thus 
predicted too high in energy. The functional also overes- 
timates the energy near equilibrium, especially with the 
empirical crpa rescaling. However, there is a definite 
improvement over both LC-wLDA and MP2. 

The addition of long-range RPA correlation to range- 
separated hybrid density functionals is a promising route 
to systematically incorporating nonlocal correlation ef- 
fects into DFT. LC-cjLDA+dRPA can be evaluated in 
0(N 4 ) time via Cholesky decomposition of A and B. 
It is robust to unstable and degenerate reference states, 
and shows promise for systems from metals, to covalent 
bonds, to van der Waals complexes. 
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